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Abstract 

A posteriori error estimates are derived in the context of two-dimensional structural elastic 
shape optimization under the compliance objective. It is known that the optimal shape fea¬ 
tures are microstructures that can be constructed using sequential lamination. The descriptive 
parameters explicitly depend on the stress. To derive error estimates the dual weighted resid¬ 
ual approach for control problems in PDE constrained optimization is employed, involving 
the elastic solution and the microstructure parameters. Rigorous estimation of interpolation 
errors ensures robustness of the estimates while local approximations are used to obtain fully 
practical error indicators. Numerical results show sharply resolved interfaces between regions 
of full and intermediate material density. 


1 Introduction 

In this paper we will address error estimation for shape optimization problems and study an 
adaptive finite element method based on an optimal microstructured elastic material composed of 
sequential laminates. A classical shape optimization problem amounts to distribute a fixed amount 
of rigid material throughout a domain in order to optimize a particular cost functional when con¬ 
fronted with surface and volume loads. However it is known that the problem is ill-posed if only 
scale invariant cost functionals are taken into account. As the formation of microstructures comes 
along with this deficiency most approaches enforce regularity by an additional regularizing cost 
functional such as the shape perimeter. Alternatively the problem can be relaxed to allow for more 
general shapes characterized by local effective material properties and densities. They stem from 
the shape structure on the microscale influencing the effective macroscopic description via homog¬ 
enization. Optimal microstructures exist but are not unique. A particular flexible microstructure 
is constructed via sequential lamination. This allows for a convenient algorithmic treatment as a 
small set of parameters characterizing the local microstructure can be computed explicitly from 
the local stress. 

Numerical computations suggest that optimal shapes are characterized by regions with substan¬ 
tially different mircostructure and sharp transitions between these regions. Thus an adaptive mesh 
strategy based on an a posteriori error control is evidently the appropriate strategy. To this end, 
the error in the achieved cost has to be controlled. 

Review of related work. Shape optimization is a classical field in PDE constrained optimization 
and has been covered extensively in the literature, see e. g. the textbooks [Ben951 rA1102] . For the 
link to homogenization theory we refer to |,IK()94l[BnMll(lD99irBDM1 . Optimal microstructures 
where first derived by Hashin in 1962 via the concentric sphere construction for hydrostatic loads 
|Has62| . The sequential lamination construction dates back to the 1980s |Tar85[ IMT851IFM86I 
ILC861IGC871 IAve87| and was later used in a practical numerical scheme for topology optimization 
in [ABF.I97] . In all these cases proofs of optimality rely on the Hashin-Shtrikman bounds on the 
attainable sets of effective elastic properties [HS63] . Related to the homogenization approach is the 
so called free material optimization, see e. g. [HKLSIO] . Here the coefficients of the elasticity are 
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the degrees of freedom to be optimized. Possibly yielding materials without a physical equivalent, 
best approximating realizations are searched in a post processing step, sometimes referred to as 
inverse homogenization. Applying the classical homogenization paradigm one asks for the optimal 
(geometric) microstructures and the effective material properties such that a macroscopic cost func¬ 
tional is optimized. In [BTlOal IBTlObj a displacement approach on the microscale is investigated 
which consists of afhne plus periodic functions on generalized periodic domains. For the numerical 
optimization of the microscopic geometry a boundary tracking method was used. The incorpora¬ 
tion of topological derivatives allowed for the nucleation of small holes and additional remeshing 
steps ensured mesh quality. For specific geometric microstructures the effective stress response of 
basic affine strains on the microscopic cell boundary are investigated in [BT12) and used to describe 
the effective macroscopic behavior. The geometry of the locally periodic microstructures are then 
optimized. The same methodology was picked up in |CGRS14] while investigating microstructured 
materials with geometrically simple perforations. The resulting shapes offered compliance values 
close to optimal ones generated by twofold sequential lamination. Moreover they indicate the need 
to resolve sharp interfaces between regions with different types of local microstructure. To be able 
to both capture microscopic effects and macroscopic material behavior while maintaining feasible 
computational complexity in PDF simulation tailored numerical methods have to be proposed. The 
heterogeneous multi-scale method (HMM) presented in |EEH03l IEE03[ IEE051IEMZ05) depicts a 
very general paradigm using independent macroscopic and microscopic schemes. In [OhlOSi [11009] 
a posteriori error estimates for elliptic homogenization problems involving a hne scale diffusion 
term were derived. Here the heterogeneous multi-scale method is reformulated as a direct finite 
element discretization of the two-scale homogenized equation in variational form. The obtained 
error indicator provides separate terms for errors in the macroscopic and microscopic discretization. 
Concerning the a posteriori control it is often not desirable to measure errors in classical function 
space norms but the error for a prescribed cost functionals. From a practical point of view quanti¬ 
ties of interest for a real composite work piece like e. g. stresses in critical sections, averaged surface 
tensions or (mollihed) pointwise stresses can be considered as cost functionals and corresponding 
a posteriori error estimates are derived |P099[ lOVOOai lOVOObi IVem04] . 

The approach we pick up in this contribution is the dual weighted residual (DWR) method intro¬ 
duced in |BR97| and applied to optimal control problems for instance in |BKR OOj . Here residual 
type error estimates get weighted using duality methods. In [BETll] a special marking strategy 
was employed allowing to show quasi-optimality of the associated adaptive finite element method 
while maintaining sound convergence. Control and state constraints were taken into account in 
|BV09[ IVW081 ILM V1 3| . Matrix valued optimal control in coefficients is studied in |KL13) . 
Shape optimization was addressed in |KV13) using a tracking type cost functional, a Helmholtz 
state equation, and a graph representation of the boundary. For this control function higher 
regularity could be shown and a priori estimates were derived. An adaptive finite element method 
for shape optimization via boundary tracking, remeshing and perimeter penalization was presented 
in [MNP V I'd! [MNPVT2]. Here a DWR type of approach is used to assess the PDE error and 
the actual geometric error is estimated differently. Furthermore the DWR method was used in 
[KDLRG141 IKall3j in the context of one shot methods for fuel ignition problems, the viscous 
Burgers equation and aerodynamic shape optimization. Optimal design in Navier Stokes flow is 
considered in [BLIJUl^ . Here the DWR method for optimal control problems leads to separate 
terms for the primal, the adjoint and the control residual. Similarly [WollOj addresses a problem 
in free material optimization based on the variable thickness sheet model. 

The paper is organized as follows. In section we will describe in condensed form the necessary 
background on linearized elasticity, shape optimization, the sequential laminates construction in 
shape optimization, and its numerical discretization. A posteriori error estimates based on the 
dual weighted residual approach are derived in sectionIn section]^ we discuss the estimation of 
weighting term using a priori regularity while section focuses on their numerical approximation. 
Details concerning the implementation are given in section and numerical results are finally 
presented in section [Tj 


2 










































2 Elastic shape optimization and microscopic lamination 

We consider an elastic object given as a simple connected domain 51 C I? c where H is a 
circumjacent working domain and we at first suppose that D is filled with elastic material (hard 
and soft) and we aim at the optimization of the distribution of the two phases. A Dirichlet boundary 
part rD C dD n 951 and a Neumann boundary Fat C dD n 951 subject to a sufficiently regular 
surface load g are both supposed to be fixed and not subject to optimization. The remaining 
boundary 951 \ (F^) U F^r) can be varied. The applied forcing causes stresses inside the material 
which maintain an inner equilibrium. It can be described by the system of partial differential 
equations of linearized elasticity. Let S denote the set of feasible stresses in accordance to the 
PDE, i.e. 


E:= {tTiZl—I div{tT} = 0 in D, an = g on Fat, u = 0 on Fd, a = C[q\£[u\] (1) 

Here u : D ^ denotes the displacement, n the outward pointing normal on Fat, £[u] = 
^ (Du + Du^) the symmetrized strain tensor, and C[q] : H —E^ the elasticity tensor. The 
elasticity tensor is a forth order tensor with the symmetry properties {C[q\)ijki = (C'[9])jifcz = 
{C[q])ijik = {C[q])kiij and the ellipticity condition C[q] ^ ^ > c^-^. Usually it can be characterized 
locally by a small parameter vector q. We consider here material which behaves isotropic and is 
described by Lame parameters A and g and the strain stress relation cr = 2g£[v\ + Atr£[u]l. 

For later purposes we take into account the weak displacement formulation of Q, i.e. M G 
the space of vector valued functions with weak derivatives in L^(51) and vanishing trace on Fa), 
solves the variational problem 


a{q\u,Lp) = l{(f) VifiG 


( 2 ) 


with the quadratic form a{q; u, (p) = C[q\ e[u] : £[(p\ dx and the linear form l{ip) = C[q\ £[u]rL- 
p da(x) = g-ip da(x). 

A classical shape optimization problem now amounts to the task of distributing a hard material 
whose elastic behavior is described by an elasticity tensor A. To this end, we consider a char¬ 
acteristic function y G L°°{D, {0,1}) of the hard phase with y dx = 0 reflecting the volume 
constraint. The elasticity tensor on D is then defined as C[q\ = X-^[q\ + ~ x)B[q] where B[q] 

is an elasticity tensor for a soft material filling the remaining space. To assess the performance 
of a given shape for a prescribed loading a suitable cost functional has to be taken into account. 
We restrict ourselves to the compliance cost as a global measure of rigidity. For a displacement 
u = u[x\ being the unique solution of © for given y it is given by 



g ■ u[x] da(x) = l{u[x]) = a{q;u[x],u[x]) ■ 


(3) 


The cost functional y] should thereby be minimized with respect to y subject to a constraint 

on the volume of the hard material phase. Ultimately one is interested in the degenerate case B — 0, 
where there is void outside of 51 = {x G U | y(a;) = 1}. This minimizing problem is ill-posed and 
on a minimizing sequence the onset of microstructures can be observed. Thus one investigates a 
relaxed formulation of this shape optimization problem. We refer to |ABF.T97llA1102| for a concise 
summary and references therein for further details. Relaxation leads to the following reformulation: 


min / min min C*\q\ ^a : a A19 A:x.. (4) 

7p,o<e<i c*[9]GGf 

Reading from left to right, for fixed stresses tr on Zl solving Q, a fixed point x G D and a fixed 
local density 0{x)^ an elasticity tensor C*[q] minimizing the elastic energy density is to be found 
among all effective homogenized elasticity tensors resulting from a mixture of materials A and B. 
Additionally the Lagrangian multiplier I takes into account the amount of total material used. 
It can be shown that it is a decreasing function of I which facilitates an adaptation to enforce 
a fixed prescribed volume at any time. For the set Gf of homogenized tensors, also referred to 
as G-closure, no algebraic characterization is known. However there exist sharp bounds on the 
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elastic energy density, called Hashin-Shtrikman bounds. Theses bounds are attained within the 
subclass of sequentially laminated materials. The construction starts on a scale ei by layering rigid 
material A and soft material B with a certain ratio and along a certain direction. By means of 
homogenization effective material properties can be computed and used in the next iteration in 
place of material B. Mechanically this takes place on a substantially larger scale e 2 ^ ei. In two 
dimensions and for the compliance cost it is known that the construction on these two iterations 
is sufficient, cf. |KL88[ fAK93b| . It is moreover possible to pass from the weak material B to void 
|AK98allXBFj^ . 

The sequential laminates microstructure is locally characterized by three parameters (cf. Fig. [^: 
the inclination of the main lamination direction a(x), the ratio of material spent in the second 
lamination stage m(x), and the overall local material density 0{x). The second lamination direction 
is always orthogonal to the first one and the second ratio parameter is 1 — m(x). 



Figure 1: Sketch of twofold sequential lamination. Numerically optimized result for an L-shaped 
domain fixed at the bottom with downwards pointing load at the center region on the right hand 
side. An adaptively refined grid is shown along with the density lamination parameter 0 plotted 
and the von Mises stress using the color coding 


All three parameters directly depend on the eigenvalues Ai, A 2 and the eigenvectors 

of the local stress tensor a and can be computed explicitly via the mapping x 1 —>■ q[u\{x) = 

{a[a\{x),m[a][x),9[a]{x))^, see |ABFJ97| : 


r 1 i. ( 

a[a\ = arctan —- , 

\tx To-)/ 


|A2(a)| 


|Ai(a)| + |A2(a)|’ 


9[a] = min < 1, 


/ 2/r + A 
4/r(^ + A) / 


(5) 


(|Ai(a)| + |A2(a)|) 


Here A and fi denote the Lame parameters of the rigid isotropic material A and , Ty'^ 
and second component of the two dimensional vector Given these parameters the 
homogenized elasticity tensor C* [g] (x) can likewise be computed explicitly: 


the first 
effective 


^mnop[^] — 7?[q;] .— t/m/o:] f/nj [o:] Qofc [ck] Qp/[o:] C ijf^l (^Ul^ ^ 


9] 

C2222[w, 9] 
Cii22[m, 9] 


4«:^(k + fj,)9{l — 9{1 — m)){l — m) 
— m)9^ + (k + A^)^(I — 9) ’ 
4K/i(K + fj,)9{l — 9m)in 
4«:^m(l — + (k + /i)^(I — 9) ’ 

4K/rA0^m(l — m) 

4K^m(l — m)9‘^ + (/c + /r)^(I — 9) 


( 6 ) 


Thereby R denotes the linear mapping which transforms the tensor C given in reference config¬ 
uration into the appropriate coordinate frame spanned by the eigenvectors of the stress tensor. 
In its definition Q are 2x2 rotation matrices and the Einstein summation convention is used. 
The bulk modulus is defined as /c = A -f /r. The tensor C is complemented using the symmetry 
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relations and filling the remaining entries with 0. This yields a singular elasticity tensor in the 
sense that a strain with just off diagonal entries, corresponding to shearing deformations, will lead 
to a vanishing stress. 

In the following we will assume that the density 9 is always bounded away from 0 by a small 
threshold. Thus we deal with the hardsoft configuration. For the ease of notation we extend the 
Neumann boundary to T^y = dD \ Tjj and prescribe homogeneous boundary conditions 5 = 0 on 
the added parts. 

For the discretization we use a finite element ansatz and assume that we are given a triangulation 
T oi D with elements T G T and denote by h the piecewise constant mesh size function. Let 
^ be the space of piecewise bi-quadratic and continuous vector valued functions on which we 
ask for a solution of the discrete weak problem a{qh\u^, ip^) = li^ph) for all G % or more 
explicitly 

[ C^lQh] e[uh] : e[iph] dx = /* g ‘ (ph da(x) ph ^ . (7) 

Jd 

Associated to Uh are stresses ■= cF[uh\{x) = C*[qh]£[uh] for which lamination parameters are 
computed via formula ([^. For brevity of notation we write qh = {ah,mh,9h) := q[crh]- Let us 
emphasize that the parameter function varies from point to point. More details on the numerical 
realization will be given in section]^ 


3 A posteriori error estimates based on the dual weighted 
residual approach 

In this section we will derive estimates of the error in the compliance objective ^ when using 
optimal sequential laminated microstructures in shape optimization. Let u be the solution of the 
continuous problem involving the elasticity tensor C* [g] with optimal lamination parameters q 
defined by ([^ and Uh be a discrete solution satisfying 0 . 

We are interested in the a posteriori control of the numerical approximation error in the cost 
functional. To this end we pick up the dual weighted residual approach [BKIIOOj . In that respect 
we will treat the rotation parameter a differently from the other two parameters m and 0 , which 
will from now on be denoted by i!) := (m, 9). 

Let us begin with a useful observation for the elastic energy density, namely that it is invariant 
with respect to the rotation of the microstructure. Intuitively this follows from the fact that the 
lamination construction yields an orthotropic material which is always aligned with the main stress 
directions. 

Lemma 3.1 (Invariance of the elastic energy density w. r. t. rotations). Let x G D be fixed and 
C*[a^'d] be an effective elasticity tensor corresponding to an optimal rank 2 sequential laminate at 
X. Let be the solution of (H). Then the local elastic energy density is invariant w. r. t. the 

rotation parameter a, i. e. 

J 

{C* [a, d] e[rt[a, dll : e\u\a^ d]]) (x) = 0 . 

da 

Proof. With (T[a,d] := C*[a,d]e[u] the elastic energy density can be rewritten as 

C*[a, d] £[u[a, d]] : e[u[a, d]] = C* ^[a, d] tT[a, d] : a\a, d]. 

Here a and d enter the effective homogenized tensor C'*[a,d] by means of Since C'*[a,d] rep¬ 
resents an optimal material the elastic energy density attains the lower Hashin-Shtrikman bound, 
cf. [A11021 Theorem 2.3.35], and we obtain 

C'*"^[a,d] cr[a,d] : (T[a,d] = A"V[a,d] : cr[a,d] -F ■ 

Here the eigenvalues Ai, A 2 of a do not depend on the rotation. In fact, the rotation parameter a 
is derived from the orientation of the eigenvectors. Thus a only appears via a in the first term of 
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the right hand side. This term also represents an elastic energy density but this time involving the 
isotropic constituent A. Because of the isotropy of A altering the alignment of ct, e. g. by rotating 
to reference configuration, leaves this term unchanged. Altogether this means that the right hand 
side does not depend on the rotation parameter. This proofs the claim. □ 


We now present the main result. 

Theorem 3.2 (Weighted a posteriori error estimate). For a given continuous solution u to with 
a microstructured material given by the optimal lamination parameter function q and a numerical 
approximation Uh the following error estimate holds for the compliance objective ([^.' 

\J[q]u] - J[qh]Uh]\ < '^r]T{uh,qh) + . 

T 


Here TZ is a higher order remainder term involving higher order derivatives and the local error 
indicators tjt of first order are defined by 

„ \_ Au) (u) (u) (u) 1 (m) (m) i (6) (6) . , 

qh) •— P'j' \ Pq'jp^q'j' ~r ^Pt 2 Pt wxtri 


Pt — 

(“) _ 
PdT — 


Pt ~ 
(S) 

Pt = 


||-div{cr?,}||o 2 ,r > 

f II5MII0.2.OT; dm 30 = % 

\ \\(Xhn — g\\o,2,dT] OTcTm ’ 

\\R[ah\C^rn[mh,Oh] e[uh[ah\] : e[uh[ah]]\\Qj^j, , 
\\R[ah]C^g[mh.,Oh] e[uh[ah]] : e[uh,[ah]]\\^ .^ .j, , 


( 2 ) 

where Uh&Vf^ 


and rfih, dh o-re chosen arbitrarily. 


4“> = 
4r’ = 
4”’ = 


1 

U-Uh 

lo, 2 ,T ’ 

1 

U-Uh 

lo,2,aT ’ 


4»> = 

e-k 

0,oo,T 


Proof. The first step is to use the Lagrangian L(q\ u,p) := J[q; u] + a(g; u,p) — l{p) for decoupled 
q, u, and p and take into account the PDE constraint, namely that the displacement u has to be 
a solution of the linearized elasticity system ([^. The first order optimality conditions are 


p[q\u){v) := £^p{q\u,p){v) = a{q;u,v) - l{v) , (8) 

p*{q;u,p)(.v) := C^u{q-,u,p){v) = J^u[q;u]{v) + a{q-,v,p) , (9) 

p'‘{q;u,p){v) := C^g{q;u,p){v) = J^q[q;u]{v) + a^q{q;u,p){v) - l,q{p){v) , (10) 


where u is a universal variable for an arbitrary test function from the corresponding space. Equation 
^ is just the definition of the elastic solution u. Equation ([^ defines the dual solution, which in 
case of compliance is p = —u. As we assumed u and Uh to be solutions of the continuous and the 
discrete problem, respectively, {q, u, —u) and (g/j, Uh, —Uh) are stationary points of the Lagrangian. 
Thus we have 


C{q]U,p) - C{qh-,Uh,Ph) = J[q]u] - J[qh-,Uh] ■ 

Estimating the error can therefore be done by considering the difference in the Lagrangian. 


Next we consider a suitable first order expansion of the Lagrangian. Here we will take special 
care of the rotation parameter a to be able to use Lemma |3.1| above. The Lagrangian thus only 
depends on u and the controls q = (a, d) because the dual solution coincides with the negative 
primal solution. Hence, in the following we will skip the dual solution as a parameter. We now 
study the dependence of the elastic solution u[a\ on the rotation parameter a. 

As a shortcut notation we write eu(s) = u[a\ — Uh[ah + sCa], Ca = ct — ah, and eg = d — 
dh for the difference between the exact and the discrete displacement and controls. Here Uh is 
given for the continuous spatially varying parameter function ah + sCa as the solution to the 
corresponding problem 0 . We now express the error in the Lagrangian as a ID integral and 
observing Uh[ah + sea] + seu(s)|^^^ = u[a] and Uh[ah + sea] + se„(s)|g^p = Uh[ah\ we obtain 

ec ■■= C{a,d-,u[a]) - C{ah,'dh-,Uh[ah]) 
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Now, we use the trapezoidal rule /(s)ds = ^(f(0) + /(I)) + 5 Jq f"(s)s(l — s) ds with f(s) = 
^£{ah + SBoc'dh + se.^,Uh[ah + sCa] + se„(s)). The term /(I) vanishes because (a,'d,M[a]) is 
assumed to be a stationary point and e = {ea,e^,u[a\ — Uh\a\f" is a feasible test direction, i.e. 
V£(a,d,u[a])(e) = 0. Thus, we obtain 


1 / Cq 

ec = -V C{ah,'dh,Uh[ah\) 

\ Uh^a[<^h]{ea) + e„(0) 


+ 7^ 


where TZ := ^ £sC{ah + seajdh + ses, Uh[o:h + sCa] + seu(s)) s (1 — s)ds . Next, we expand the 
lower order term and achieve 


1 / Cq 

li)V,fi{oih\') I 

\ Uh^a[<^h]{ea) + e„( 0 ) 

= ^C^u{a:h,'dh;uh[ah]){u[a] - Uh[ah]) + ^C^^{ah,dh]Uh[ah\){'d - -dh) 

+ ^^,a{ah,dh':Uh[ah]){a - ah) + ^C^u{oih,^h;uh[ah]){uh^a[ah]ia - ah)) 

= ^^,uiah,'&h;uh[ah]){u[a] - Uh[ah]) + ^C^i}{ah,^h;uh[ah])id - dh) 

1 d 

By definition we have 


d 

da 


C{a,'dh;uh[a]) 


ct—ah 


d 

da 


'D 


C* [a,'dh]s[uh[a]] : e[uh[a]] dx 


and thus the term vanishes by virtue of Lemma [3.1[ Finally we consider arbitrary discrete test func¬ 
tions Uh G and dh resulting from an arbitrary stress by means of formula ([^ and use the fact 
that {ah,'&h,Uh[ah]) is a discrete stationary point of the Lagrangian and thus VC{ah, 'dh, Uh[cth]) 
vanishes in directions (0,'dft, — dh,Uh[oi-h\ — Uh[oih\)- This gives 

ec = ]^C,u{oih,'dh-,Uh[ah\){u[a\ - Uh[ah\) + ]^C^^{ah,'dh]Uh[ah\){'d -^h) d-TZ- 

Furthermore, we consider the weak formulation ([^ and the compliance cost functional (|^. For 
the primal residual we get: 


d^Ao^h,'dh]Uh[ah])iu[a] - Uh[o^h]) 

= -2/ C* [ah,'dh]e[uh[ah]] ■ e[u[a\ - Uh[ah]] dx + 2 g ■ {u[a\ - Uh[ah]) 
J D JTn 

/ div {C* [ah, id h]£[uh[ah]]} ■ {,u[a]-Uh[ah]) d:x. 


= 2 


IdTi 


C* [ah, d/,] e[uh[ah]]n ■ {u[a] - Uh[ah]) da(x) 


g ■ {u[a] - Uh[ah]) da(x)^ 


dTj nr j\j- 

This leads to the postulated residual terms with 


= ||div{CT[u/,K]]}||o 2 ,T> 

>) _ / HWiMo^h\]-n]\\o, 2 ,dT 
I \W[uh[ah]]-n-g\\Q^ 2 ,aT 


dTj n az? = 0 
dTj c Fat 


uP = ||u[a] -u/*K]||o 2 ,t> 




For the control residual we get using Q 
^,^{ah,^h;uh[oih]){^ - -^h) 

= V / [ah,-dh] (t^ - '^h)£b^h[ah]] ■ £[uh[ah]] dx 
= ^ / R[ah] Cm [mh,0h] {m - rhh) e[uh[ah]] ■ e[uh[ah]] 

j 

+R[ah] Cg [mh,6h] {9 - Oh) e[uh[ah]] : £[uh[ah]] dx 

Concerning the last expression we know that the controls ruh. Oh being the laminate parameters 
are bounded. Differentiation of the entries of the effective tensor Q w. r.t. m and 0 again yields 
bounded expressions within the attainable range of the parameters. The tensors C\m and C^g are 
therefore bounded and the whole integrand is bounded in Li. We finally obtain the desired residual 
estimate with 


= \\R[ah]C 

,m [mh,9h] £[uh[ah]] : £[uh[ah]]\\Q j^ j, , 

piP = \\R[ah]C^e[mh,Oh] £[uh[ah]] : e[uh[ah]]\\f^^j^ ,j, , 




Wt = - ^h\\o,oo,T 

O-Oh 


, ,(e) _ 

(jjrp - 


0,oo,T 


This proofs the claim. 


□ 


4 Estimation of weights using a priori regularity 

For the error estimates of Theorem |3.2| the residual terms px can be computed straightforwardly 
from the discrete solution Uh- The weights however, depend on the exact (unknown) solution 
u. Using a priori known regularity those terms can be estimated. For the primal solution classical 
interpolation estimates exist. For the laminate parameters it is possible to estimate the error by 
exploiting the explicit relation to the primal solution ([^. The following estimates especially ensure 
robustness of the error estimates for a sufficiently smooth continuous solution u: 

Corollary 1 (A priori estimates of the weights). Under the assumption that u G the 

following estimates for the weighting terms hold: 


I — 

UJrp - 

, ,(“) _ 
LOqt — 

, ,(™) _ 
UJrp - 


UJrp - 


11 '*^ '^^llo,2,T 
IIm ~ '^/i|lo,2,aT 




9-k 


0,oo,T 


< h{T) \u\i^2,ujt ) 

< h{T)i \u\i^2,ujaT : 
<h{T)\\nh{T)\ |u| 2 ,oo.r, 
<h{T)\\nh{T)\ |M| 2 ,oo.r. 


( 11 ) 


Here, h{T) is the edge length of the square element T and ^ denotes the smaller or equal inequality 
up to a constant factor depending solely on the triangulation. 

Proof. For the primal error weights we choose Uh = ^h where denotes the piecewise bi¬ 
quadratic Lagrangian interpolation. Then standard interpolation estimates |Cia78[ Theorem 3.1.5] 
give 


||m — T^^^u||o,2,T ^ h{T)\u\ip^T , 

||m — dJh'’ u\\Q^2,dT ^ h{T)^\u\i^2,dT ■ 

Next, we estimate the error with respect to the control parameter m. The estimate for 0 is 

(2) 

completely analogous. We choose rhh = m[C*£[If 'u]] =: in{Xi{ah), X 2 {<^h)), i-e. rhh is obtained 
by the pointwise evaluation of the interpolated solution u, the computation of the corresponding 
stress (T/j = C'*e[I^^^it]] and its eigenvalues and finally inserting them into formulae Then, 
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sup^gT^ |m(Ai((T), A 2 (cr))— m(Ai(cr^), A 2 (cr^))| is attained for some x G T with stresses cr = cr{x), 
= c^hix). From the Lipschitz continuity of the function m we deduce 


m 


(Ai(cr), A 2 (fT)) - m(Ai((T,i), A 2 (cr?i))| < (|Ai(cr) - Xi{ah)\^ + |A 2 (cr) - X 2 iah)\ 


2 \ 2 


Furthermore, the eigenvalues are Lipschitz continuous functions of the underlying matrices. In 
particular the Wielandt-Hoffmann inequality |HW53j for p = 2 yields 

(|Ai(ct) - Ai(d-/i)p + |A2(d-) - A2(d-;i)P) ^ < \\a - anWr ■ 

The elasticity tensor C* represents a bounded linear operator, thus 




0,00,T 


Cie[u\ - e[xj^K]) 


< 


0,oo,T 


Mu]- 4^h^u]\\o,oc,T 


with 111(7111 representing the associated operator norm of the tensor C. Passing from the sym¬ 
metrized strain tensor to the full Jacobian yields 

||e[u] - < |u • 

Finally, by standard L°° estimates |Cia78[ Theorem 3.3.7], we achieve 

\u - Xf ^uli.oo.T < h{T)\lnh{T)\ |u| 2 ,oo,t • 


□ 

Remark 1. The presented argument would not work for a weighting term like ||o — 5/i||g ^ 
involving the rotation parameter as it depends on the eigenvectors of the local stress and the 
eigenvectors do not depend continuously on the matrix. In fact this motivates the elimination of 
the term involving the rotation parameter in Theorem |3.2[ 


5 Numerical treatment 

For a practical numerical scheme for the weighting terms ojt that were a priori estimated in section]^ 
we ask for an approximation based on the numerical solution Uh- One approach would be to exploit 
higher regularity of the solution by computing a second discrete solution using a higher order 
scheme. However this is usually computationally demanding. To keep the cost for the evaluation 
of the error estimates low |BR97| suggests to merely interpolate the discrete solution to a higher 
order polynomial. As we use bi-quadratic hnite elements to overcome checkerboard instabilities, 
cf. |JB98| , we therefore construct polynomials of degree four on patches of four neighboring elements 
from the local degrees of freedom of the discrete elastic displacement Uh- 

In the discrete probl em ([t]) the coefficient functions q and thus C* vary from point to point. How¬ 
ever it is known, see |Cia781 Theorem 4.1.6], that the convergence order of the finite element scheme 
is unaltered if the chosen quadrature order is high enough, which is the case in our computations. 
To this end the parameters a, m can easily be evaluated at the quadrature points. The same holds 
true for the density parameter as well. However it needs to be piecewise constant on each element 
to avoid checkerboard instabilities, see |JB98| . This is ensured by an element wise averaging of all 
values obtained on quadrature points. In detail we proceed as follows. In case of the ratio param¬ 
eter m (and the rotation parameter a which however was eliminated from the error estimate) the 
interpolated discrete deformation Uh can be used again to locally obtain parameters a and m for 
fixed d. For the density parameter the piecewise constant approximation can be used to construct 
a bilinear profile on a patch of four neighboring elements. Let us summarize the approximation of 
the weights in the following remark: 

Remark 2 (Approximation of local weights). To approximate the weights for the primal and the 
control error indicator on every element X and to derive a suitable error indicator we consider a 
patch S of the four neighboring elements including X and proceed as follows: 
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- For the primal error indicator 5x5 degrees of freedom on each element of the patch are 
used to construct a bi-quartic polynomial on the patch. This interpolation is then 

compared to the original approximation leading to 


(jJrj 


^(4) 

Uh 'Uh 


o,2.r 


UJ 


{u) 

dT 


^(4) 

Uh -Ih 


0,2,ST 


- For the control error indicator the ratio parameter m is computed from the interpolated 
elastic solution above and compared to the original value; for the density parameter 9 four 
piecewise constant values are used to construct a bilinear profile on the patch which is then 
compared to the original values. In explicit 


,(”i) 


m[uh] - 


,(») 


0,oo,T 


Oh-T^h^dh 


0,oo,T 


6 Implementation 

For our computations we use an adaptive regular quadrilateral mesh provided by the QuocMesh 
library!^ Refinement is done via uniform refinement of cells and handling of constrained hanging 
nodes. To obtain the optimal sequential lamination microstructure we reimplemented the alternat¬ 
ing algorithm suggested in |ABFJ97j . The effective tensors are regularized by setting C 44 = 10“^. 
To ensure coercivity we restrict m to the range [e, 1 — e] and 6 to [e, 1] with e = 10“^. The 
algorithm is terminated once the change in the compliance cost value is below 10“^. As already 
mentioned we use bi-quadratic elements to overcome checkerboard instabilities. For the numerical 
integration it turned out to be sufficient to use a Gauss quadrature rule of order 5. The lamination 
parameters for the interpolated solution cannot be directly computed from the pointwise stress 
C* £\I^^Uh] as it in turn depends on the unknown parameters. Here we use a Newton method to 
solve the equation 

C[m{\i), 9] u„] - Q{a) Q{a)^ = 0 

for the rotation parameter a and the eigenvalues Xi of the stress (from which m can be computed). 
The density 9 is fixed to the original value found at the current point. 

For the adaptive scheme we use a Dorfler marking strategy [Dbr96) with a fraction of 40% of the 
total error to be associated with the set of elements to be refined. 


7 Numerical results 

As the first example we consider a carrier plate under shearing computed on the unit square., 
cf. Fig. The volume constraint was set to 33%. First we ran the alternating algorithm until 
convergence on uniform meshes of level I G {2,..., 10}. The obtained values can be used to 
extrapolate the asymptotic value of the compliance cost by assuming the following dependence on 
the grid width: 

J[Uh] = J*+chP. 

Using a least squares fit we found the parameters J* = 1.8399, c = 1.7645, and p = 1.0484. J* is 
used to replace the exact cost value when assessing error reduction in the following. 

We ran the discussed adaptive scheme for the carrier plate scenario, see Figures and 

Next, we take into account a eantilever scenario. Here we consider the domain D = [0, 2] x [0,1] 

with Dirichlet boundary condition on the left hand side and a downward pointing force in the 

middle on the right hand side. The volume constraint is set to 50%, cf. Figure 

As a third scenario we investigate a bridge within the domain D = [0, 2] x [0,1] in Figure We 

prescribe roller boundary conditions on short strips at the lower boundary in the left and right 

^http://numod.ins.uni-bonn.de/software/quocmesh 
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Figure 2: Adaptive meshes after 4, 7, 8, 10, 13, 16, 19, and 20 refinement steps obtained for the 
carrier plate scenario. Color plots show the error indicator after 7 and 19 steps, leading to the 
subsequent grid. 



Figure 3: Difference to extrapolated value J* plotted over the number of elements, once for a 
uniform refinement (green) and the adaptive refinement strategy (red). 


corners, i. e. the displacement in j/-direction is enforced to be 0, compare for instance |A1102) . 
Additionally the node in the lower left corner is subject to full homogeneous Dirichlet boundary 
conditions to ensure a unique solution. The lower boundary in between is loaded in vertical 
direction. The volume constraint is set to 33%. 


8 Conclusion 

In this paper we have applied the dual weighted residual approach to an elastic multi scale shape 
optimization problem using sequential laminates as an underlying optimal microstructure. The 
microstructure enters the PDF via effective material properties whose explicit formulae allow to 
compute sensitivities of the cost functional with respect to the describing parameters of the mi¬ 
crostructure. The derived goal-oriented error estimate enables to appropriately steer the adaptive 
refinement leading to a substantial reduction of the number of degrees of freedom compared to uni¬ 
form meshes and the same accuracy. In particular the goal oriented error estimation outperforms 
a residual type estimator solely for the elastic displacement problem. The main ingredient is the 
error term incorporating the sensitivity of the elastic energy density w. r. t. to the local material 
density. This is in accordance with the observation that the rigidity of a shape is mostly improved 
by an optimal distribution of the available material. The adaptive scheme leads to sharply resolved 
interfaces, separating regions with almost no material from regions with bulk material or an inter¬ 
mediate material density. 
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Figure 4: At different refinement steps of the adaptive method applied to the cantilever example 
the current mesh (top) and the corresponding density 9 (bottom) are displayed (from left to right: 
step 4, 8, 12, 16, 20, and 24). 
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